25class PseudoDifferentialOperator:
  26    """
  27    Pseudo-differential operator with dynamic symbol evaluation on spatial grids.
  28    Supports both 1D and 2D operators, and can be defined explicitly (symbol mode)
  29    or extracted automatically from symbolic equations (auto mode).
  30
  31    Parameters
  32    ----------
  33    expr : sympy expression
  34        Symbolic expression representing the pseudo-differential symbol.
  35    vars_x : list of sympy symbols
  36        Spatial variables (e.g., [x] for 1D, [x, y] for 2D).
  37    var_u : sympy function, optional
  38        Function u(x, t) used in auto mode to extract the operator symbol.
  39    mode : str, {'symbol', 'auto'}
  40        - 'symbol': directly uses expr as the operator symbol.
  41        - 'auto': computes the symbol automatically by applying expr to exp(i x ξ).
  42
  43    Attributes
  44    ----------
  45    dim : int
  46        Spatial dimension (1 or 2).
  47    fft, ifft : callable
  48        Fast Fourier transform and inverse (scipy.fft or scipy.fft2).
  49    p_func : callable
  50        Evaluated symbol function ready for numerical use.
  51
  52    Notes
  53    -----
  54    - In 'symbol' mode, `expr` should be expressed in terms of spatial variables and frequency variables (ξ, η).
  55    - In 'auto' mode, the symbol is derived by applying the differential expression to a complex exponential.
  56    - Frequency variables are internally named 'xi' and 'eta' for consistency.
  57    - Uses numpy for numerical evaluation and scipy.fft for FFT operations.
  58
  59    Examples
  60    --------
  61    >>> # Example 1: 1D Laplacian operator (symbol mode)
  62    >>> from sympy import symbols
  63    >>> x, xi = symbols('x xi', real=True)
  64    >>> op = PseudoDifferentialOperator(expr=xi**2, vars_x=[x], mode='symbol')
  65
  66    >>> # Example 2: 1D transport operator (auto mode)
  67    >>> from sympy import Function
  68    >>> u = Function('u')
  69    >>> expr = u(x).diff(x)
  70    >>> op = PseudoDifferentialOperator(expr=expr, vars_x=[x], var_u=u(x), mode='auto')
  71    """
  72
  73    def __init__(self, expr, vars_x, var_u=None, mode='symbol'):
  74        self.dim = len(vars_x)
  75        self.mode = mode
  76        self.symbol_cached = None
  77        self.expr = expr
  78        self.vars_x = vars_x
  79
  80        if self.dim == 1:
  81            x, = vars_x
  82            xi_internal = symbols('xi', real=True)
  83            expr = expr.subs(symbols('xi', real=True), xi_internal)
  84            self.fft = partial(fft, workers=FFT_WORKERS)
  85            self.ifft = partial(ifft, workers=FFT_WORKERS)
  86
  87            if mode == 'symbol':
  88                self.p_func = lambdify((x, xi_internal), expr, 'numpy')
  89                self.symbol = expr
  90            elif mode == 'auto':
  91                if var_u is None:
  92                    raise ValueError("var_u must be provided in mode='auto'")
  93                exp_i = exp(I * x * xi_internal)
  94                P_ei = expr.subs(var_u, exp_i)
  95                symbol = simplify(P_ei / exp_i)
  96                symbol = expand(symbol)
  97                self.symbol = symbol
  98                self.p_func = lambdify((x, xi_internal), symbol, 'numpy')
  99            else:
 100                raise ValueError("mode must be 'auto' or 'symbol'")
 101
 102        elif self.dim == 2:
 103            x, y = vars_x
 104            xi_internal, eta_internal = symbols('xi eta', real=True)
 105            expr = expr.subs(symbols('xi', real=True), xi_internal)
 106            expr = expr.subs(symbols('eta', real=True), eta_internal)
 107            self.fft = partial(fft2, workers=FFT_WORKERS)
 108            self.ifft = partial(ifft2, workers=FFT_WORKERS)
 109
 110            if mode == 'symbol':
 111                self.symbol = expr
 112                self.p_func = lambdify((x, y, xi_internal, eta_internal), expr, 'numpy')
 113            elif mode == 'auto':
 114                if var_u is None:
 115                    raise ValueError("var_u must be provided in mode='auto'")
 116                exp_i = exp(I * (x * xi_internal + y * eta_internal))
 117                P_ei = expr.subs(var_u, exp_i)
 118                symbol = simplify(P_ei / exp_i)
 119                symbol = expand(symbol)
 120                self.symbol = symbol
 121                self.p_func = lambdify((x, y, xi_internal, eta_internal), symbol, 'numpy')
 122            else:
 123                raise ValueError("mode must be 'auto' or 'symbol'")
 124
 125        else:
 126            raise NotImplementedError("Only 1D and 2D supported")
 127
 128        if mode == 'auto':
 129            print("\nsymbol = ")
 130            pprint(self.symbol, num_columns=NUM_COLS)
 131        
 132    def evaluate(self, X, Y, KX, KY, cache=True):
 133        """
 134        Evaluate the pseudo-differential operator's symbol on a grid of spatial and frequency coordinates.
 135
 136        The method dynamically selects between 1D and 2D evaluation based on the spatial dimension.
 137        If caching is enabled and a cached symbol exists, it returns the cached result to avoid recomputation.
 138
 139        Parameters
 140        ----------
 141        X, Y : ndarray
 142            Spatial grid coordinates. In 1D, Y is ignored.
 143        KX, KY : ndarray
 144            Frequency grid coordinates. In 1D, KY is ignored.
 145        cache : bool, default=True
 146            If True, stores the computed symbol for reuse in subsequent calls to avoid redundant computation.
 147
 148        Returns
 149        -------
 150        ndarray
 151            Evaluated symbol values over the input grid. Shape matches the input spatial/frequency grids.
 152
 153        Raises
 154        ------
 155        NotImplementedError
 156            If the spatial dimension is not 1D or 2D.
 157        """
 158        if cache and self.symbol_cached is not None:
 159            return self.symbol_cached
 160
 161        if self.dim == 1:
 162            symbol = self.p_func(X, KX)
 163        elif self.dim == 2:
 164            symbol = self.p_func(X, Y, KX, KY)
 165
 166        if cache:
 167            self.symbol_cached = symbol
 168
 169        return symbol
 170
 171    def clear_cache(self):
 172        """
 173        Clear cached symbol evaluations.
 174        """        
 175        self.symbol_cached = None
 176
 177    def apply(self, u, x_grid, kx, boundary_condition='periodic', 
 178              y_grid=None, ky=None, dealiasing_mask=None,
 179              freq_window='gaussian', clamp=1e6, space_window=False):
 180        """
 181        Apply the pseudo-differential operator to the input field u.
 182    
 183        This method dispatches the application of the pseudo-differential operator based on:
 184        
 185        - Whether the symbol is spatially dependent (x/y)
 186        - The boundary condition in use (periodic or dirichlet)
 187    
 188        Supported operations:
 189        
 190        - Constant-coefficient symbols: applied via Fourier multiplication.
 191        - Spatially varying symbols: applied via Kohn–Nirenberg quantization.
 192        - Dirichlet boundary conditions: handled with non-periodic convolution-like quantization.
 193    
 194        Dispatch Logic:\n
 195        if not self.is_spatial: u ↦ Op(p)(D) ⋅ u = 𝓕⁻¹[ p(ξ) ⋅ 𝓕(u) ]\n
 196        elif periodic: u ↦ Op(p)(x,D) ⋅ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ based of FFT (quicker)\n
 197        elif dirichlet: u ↦ Op(p)(x,D) ⋅ u ≈ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ (slower)\n
 198        
 199        Parameters
 200        ----------
 201        u : ndarray
 202            Function to which the operator is applied
 203        x_grid : ndarray
 204            Spatial grid in x direction
 205        kx : ndarray
 206            Frequency grid in x direction
 207        boundary_condition : str
 208            'periodic' or 'dirichlet'
 209        y_grid : ndarray, optional
 210            Spatial grid in y direction (for 2D)
 211        ky : ndarray, optional
 212            Frequency grid in y direction (for 2D)
 213        dealiasing_mask : ndarray, optional
 214            Dealiasing mask
 215        freq_window : str
 216            Frequency windowing ('gaussian' or 'hann')
 217        clamp : float
 218            Clamp symbol values to [-clamp, clamp]
 219        space_window : bool
 220            Apply spatial windowing
 221            
 222        Returns
 223        -------
 224        ndarray
 225            Result of applying the operator
 226        """
 227        # Check if symbol depends on spatial variables
 228        is_spatial = self._is_spatial_dependent()
 229        
 230        # Case 1: Constant symbol with periodic BC (fast path)
 231        if not is_spatial and boundary_condition == 'periodic':
 232            return self._apply_constant_fft(u, x_grid, kx, y_grid, ky, dealiasing_mask)
 233        
 234        # Case 2: Spatial symbol with periodic BC
 235        elif boundary_condition == 'periodic':
 236            symbol_func = self._get_symbol_func()
 237            return kohn_nirenberg_fft(
 238                u_vals=u,
 239                symbol_func=symbol_func,
 240                x_grid=x_grid,
 241                kx=kx,
 242                fft_func=self.fft,
 243                ifft_func=self.ifft,
 244                dim=self.dim,
 245                y_grid=y_grid,
 246                ky=ky,
 247                freq_window=freq_window,
 248                clamp=clamp,
 249                space_window=space_window
 250            )
 251        
 252        # Case 3: Dirichlet BC (non-periodic)
 253        elif boundary_condition == 'dirichlet':
 254            symbol_func = self._get_symbol_func()
 255            
 256            if self.dim == 1:
 257                return kohn_nirenberg_nonperiodic(
 258                    u_vals=u,
 259                    x_grid=x_grid,
 260                    xi_grid=kx,
 261                    symbol_func=symbol_func,
 262                    freq_window=freq_window,
 263                    clamp=clamp,
 264                    space_window=space_window
 265                )
 266            elif self.dim == 2:
 267                return kohn_nirenberg_nonperiodic(
 268                    u_vals=u,
 269                    x_grid=(x_grid, y_grid),
 270                    xi_grid=(kx, ky),
 271                    symbol_func=symbol_func,
 272                    freq_window=freq_window,
 273                    clamp=clamp,
 274                    space_window=space_window
 275                )
 276        
 277        else:
 278            raise ValueError(f"Invalid boundary condition '{boundary_condition}'")
 279    
 280    def _is_spatial_dependent(self):
 281        """
 282        Check if the symbol depends on spatial variables.
 283        
 284        Returns
 285        -------
 286        bool
 287            True if symbol depends on x (or x, y)
 288        """
 289        if self.dim == 1:
 290            return self.symbol.has(self.vars_x[0])
 291        elif self.dim == 2:
 292            x, y = self.vars_x
 293            return self.symbol.has(x) or self.symbol.has(y)
 294        else:
 295            return False
 296    
 297    def _get_symbol_func(self):
 298        """
 299        Get a lambdified version of the symbol.
 300        
 301        Returns
 302        -------
 303        callable
 304            Lambdified symbol function
 305        """
 306        if self.dim == 1:
 307            x = self.vars_x[0]
 308            xi = symbols('xi', real=True)
 309            return lambdify((x, xi), self.symbol, 'numpy')
 310        elif self.dim == 2:
 311            x, y = self.vars_x
 312            xi, eta = symbols('xi eta', real=True)
 313            return lambdify((x, y, xi, eta), self.symbol, 'numpy')
 314        else:
 315            raise NotImplementedError("Only 1D and 2D supported")
 316    
 317    def _apply_constant_fft(self, u, x_grid, kx, y_grid, ky, dealiasing_mask):
 318        """
 319        Apply a constant-coefficient pseudo-differential operator in Fourier space.
 320
 321        This method assumes the symbol is diagonal in the Fourier basis and acts as a 
 322        multiplication operator. It performs the operation:
 323        
 324            (ψu)(x) = 𝓕⁻¹[ -σ(k) · 𝓕[u](k) ]
 325
 326        where:
 327        - σ(k) is the combined pseudo-differential operator symbol
 328        - 𝓕 denotes the forward Fourier transform
 329        - 𝓕⁻¹ denotes the inverse Fourier transform
 330
 331        The dealiasing mask is applied before returning to physical space.
 332        
 333        Parameters
 334        ----------
 335        u : ndarray
 336            Input function
 337        x_grid : ndarray
 338            Spatial grid (x)
 339        kx : ndarray
 340            Frequency grid (x)
 341        y_grid : ndarray, optional
 342            Spatial grid (y, for 2D)
 343        ky : ndarray, optional
 344            Frequency grid (y, for 2D)
 345        dealiasing_mask : ndarray, optional
 346            Dealiasing mask
 347            
 348        Returns
 349        -------
 350        ndarray
 351            Result
 352        """
 353        u_hat = self.fft(u)
 354        
 355        # Evaluate symbol at grid points
 356        if self.dim == 1:
 357            X_dummy = np.zeros_like(kx)
 358            symbol_vals = self.p_func(X_dummy, kx)
 359        elif self.dim == 2:
 360            KX, KY = np.meshgrid(kx, ky, indexing='ij')
 361            X_dummy = np.zeros_like(KX)
 362            Y_dummy = np.zeros_like(KY)
 363            symbol_vals = self.p_func(X_dummy, Y_dummy, KX, KY)
 364        else:
 365            raise ValueError("Only 1D and 2D supported")
 366        
 367        # Apply symbol
 368        u_hat *= symbol_vals
 369        
 370        # Apply dealiasing
 371        if dealiasing_mask is not None:
 372            u_hat *= dealiasing_mask
 373        
 374        return self.ifft(u_hat)
 375
 376    def principal_symbol(self, order=1):
 377        """
 378        Compute the leading homogeneous component of the pseudo-differential symbol.
 379
 380        This method extracts the principal part of the symbol, which is the dominant 
 381        term under high-frequency asymptotics (|ξ| → ∞). The expansion is performed 
 382        in polar coordinates for 2D symbols to maintain rotational symmetry, then 
 383        converted back to Cartesian form.
 384
 385        Parameters
 386        ----------
 387        order : int
 388            Order of the asymptotic expansion in powers of 1/ρ, where ρ = |ξ| in 1D 
 389            or ρ = sqrt(ξ² + η²) in 2D. Only the leading-order term is returned.
 390
 391        Returns
 392        -------
 393        sympy.Expr
 394            The principal symbol component, homogeneous of degree `m - order`, where 
 395            `m` is the original symbol's order.
 396
 397        Notes:
 398        - In 1D, uses direct series expansion in ξ.
 399        - In 2D, expands in radial variable ρ while preserving angular dependence.
 400        - Useful for microlocal analysis and constructing parametrices.
 401        """
 402
 403        p = self.symbol
 404        if self.dim == 1:
 405            xi = symbols('xi', real=True, positive=True)
 406            return simplify(series(p, xi, oo, n=order).removeO())
 407        elif self.dim == 2:
 408            xi, eta = symbols('xi eta', real=True, positive=True)
 409            # Homogeneous radial expansion: we set (ξ, η) = ρ (cosθ, sinθ)
 410            rho, theta = symbols('rho theta', real=True, positive=True)
 411            p_rho = p.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 412            expansion = series(p_rho, rho, oo, n=order).removeO()
 413            # Revert back to (ξ, η)
 414            expansion_cart = expansion.subs({rho: sqrt(xi**2 + eta**2),
 415                                             cos(theta): xi / sqrt(xi**2 + eta**2),
 416                                             sin(theta): eta / sqrt(xi**2 + eta**2)})
 417            return simplify(powdenest(expansion_cart, force=True))
 418                       
 419    def is_homogeneous(self, tol=1e-10):
 420        """
 421        Check whether the symbol is homogeneous in the frequency variables.
 422    
 423        Returns
 424        -------
 425        (bool, Rational or float or None)
 426            Tuple (is_homogeneous, degree) where:
 427            - is_homogeneous: True if the symbol satisfies p(λξ, λη) = λ^m * p(ξ, η)
 428            - degree: the detected degree m if homogeneous, or None
 429        """
 430        from sympy import symbols, simplify, expand, Eq
 431        from sympy.abc import l
 432    
 433        if self.dim == 1:
 434            xi = symbols('xi', real=True, positive=True)
 435            l = symbols('l', real=True, positive=True)
 436            p = self.symbol
 437            p_scaled = p.subs(xi, l * xi)
 438            ratio = simplify(p_scaled / p)
 439            if ratio.has(xi):
 440                return False, None
 441            try:
 442                deg = simplify(ratio).as_base_exp()[1]
 443                return True, deg
 444            except Exception:
 445                return False, None
 446    
 447        elif self.dim == 2:
 448            xi, eta = symbols('xi eta', real=True, positive=True)
 449            l = symbols('l', real=True, positive=True)
 450            p = self.symbol
 451            p_scaled = p.subs({xi: l * xi, eta: l * eta})
 452            ratio = simplify(p_scaled / p)
 453            # If ratio == l**m with no (xi, eta) left, it's homogeneous
 454            if ratio.has(xi, eta):
 455                return False, None
 456            try:
 457                base, exp = ratio.as_base_exp()
 458                if base == l:
 459                    return True, exp
 460            except Exception:
 461                pass
 462            return False, None
 463
 464    def symbol_order(self, max_order=10, tol=1e-3):
 465        """
 466        Estimate the homogeneity order of the pseudo-differential symbol in high-frequency asymptotics.
 467    
 468        This method attempts to determine the leading-order behavior of the symbol p(x, ξ) or p(x, y, ξ, η)
 469        as |ξ| → ∞ (in 1D) or |(ξ, η)| → ∞ (in 2D). The returned value represents the asymptotic growth or decay rate,
 470        which is essential for understanding the regularity and mapping properties of the corresponding operator.
 471    
 472        The function uses symbolic preprocessing to ensure proper factorization of frequency variables,
 473        especially in sqrt and power expressions, to avoid erroneous order detection (e.g., due to hidden scaling).
 474    
 475        Parameters
 476        ----------
 477        max_order : int, optional
 478            Maximum number of terms to consider in the series expansion. Default is 10.
 479        tol : float, optional
 480            Tolerance threshold for evaluating the coefficient magnitude. If the coefficient is too small,
 481            the detected order may be discarded. Default is 1e-3.
 482    
 483        Returns
 484        -------
 485        float or None
 486            - If the symbol is homogeneous, returns its exact homogeneity degree as a float.
 487            - Otherwise, estimates the dominant asymptotic order from leading terms in the expansion.
 488            - Returns None if no valid order could be determined.
 489    
 490        Notes
 491        -----
 492        - In 1D:
 493            Two strategies are used:
 494                1. Expand directly in xi at infinity.
 495                2. Substitute xi = 1/z and expand around z = 0.
 496    
 497        - In 2D:
 498            - Transform the symbol into polar coordinates: (xi, eta) = rho*(cos(theta), sin(theta)).
 499            - Expand in rho at infinity, then extract the leading term's power.
 500            - An alternative substitution using 1/z is also tried if the first method fails.
 501    
 502        - Preprocessing steps:
 503            - Sqrt expressions involving frequencies are rewritten to isolate the leading variable.
 504            - Power expressions are factored explicitly to ensure correct symbolic scaling.
 505    
 506        - If the symbol is not homogeneous, a warning is issued, and the result should be interpreted with care.
 507        
 508        - For non-homogeneous symbols, only the principal asymptotic term is considered.
 509    
 510        Raises
 511        ------
 512        NotImplementedError
 513            If the spatial dimension is neither 1 nor 2.
 514        """
 515        from sympy import (
 516            symbols, series, simplify, sqrt, cos, sin, oo, powdenest, radsimp,
 517            expand, expand_power_base
 518        )
 519    
 520        def preprocess_sqrt(expr, freq):
 521            return expr.replace(
 522                lambda e: e.func == sqrt and freq in e.free_symbols,
 523                lambda e: freq * sqrt(1 + (e.args[0] - freq**2) / freq**2)
 524            )
 525    
 526        def preprocess_power(expr, freq):
 527            return expr.replace(
 528                lambda e: e.is_Pow and freq in e.free_symbols,
 529                lambda e: freq**e.exp * (1 + e.base / freq**e.base.as_powers_dict().get(freq, 0))**e.exp
 530            )
 531    
 532        def validate_order(power, coeff, vars_x, tol):
 533            if power is None:
 534                return None
 535            if any(v in coeff.free_symbols for v in vars_x):
 536                print("⚠️ Coefficient depends on spatial variables; ignoring")
 537                return None
 538            try:
 539                coeff_val = abs(float(coeff.evalf()))
 540                if coeff_val < tol:
 541                    print(f"⚠️ Coefficient too small ({coeff_val:.2e} < {tol})")
 542                    return None
 543            except Exception as e:
 544                print(f"⚠️ Coefficient evaluation failed: {e}")
 545                return None
 546            return int(power) if power == int(power) else float(power)
 547    
 548        # Homogeneity check
 549        is_homog, degree = self.is_homogeneous()
 550        if is_homog:
 551            return float(degree)
 552        else:
 553            print("⚠️ The symbol is not homogeneous. The asymptotic order is not well defined.")
 554    
 555        if self.dim == 1:
 556            x = self.vars_x[0]
 557            xi = symbols('xi', real=True, positive=True)
 558    
 559            try:
 560                print("1D symbol_order - method 1")
 561                expr = preprocess_sqrt(self.symbol, xi)
 562                s = series(expr, xi, oo, n=max_order).removeO()
 563                lead = simplify(powdenest(s.as_leading_term(xi), force=True))
 564                power = lead.as_powers_dict().get(xi, None)
 565                coeff = lead / xi**power if power is not None else 0
 566                print("lead =", lead)
 567                print("power =", power)
 568                print("coeff =", coeff)
 569                order = validate_order(power, coeff, [x], tol)
 570                if order is not None:
 571                    return order
 572            except Exception:
 573                pass
 574    
 575            try:
 576                print("1D symbol_order - method 2")
 577                z = symbols('z', real=True, positive=True)
 578                expr_z = preprocess_sqrt(self.symbol.subs(xi, 1/z), 1/z)
 579                s = series(expr_z, z, 0, n=max_order).removeO()
 580                lead = simplify(powdenest(s.as_leading_term(z), force=True))
 581                power = lead.as_powers_dict().get(z, None)
 582                coeff = lead / z**power if power is not None else 0
 583                print("lead =", lead)
 584                print("power =", power)
 585                print("coeff =", coeff)
 586                order = validate_order(power, coeff, [x], tol)
 587                if order is not None:
 588                    return -order
 589            except Exception as e:
 590                print(f"⚠️ fallback z failed: {e}")
 591            return None
 592    
 593        elif self.dim == 2:
 594            x, y = self.vars_x
 595            xi, eta = symbols('xi eta', real=True, positive=True)
 596            rho, theta = symbols('rho theta', real=True, positive=True)
 597    
 598            try:
 599                print("2D symbol_order - method 1")
 600                p_rho = self.symbol.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 601                p_rho = preprocess_power(preprocess_sqrt(p_rho, rho), rho)
 602                s = series(simplify(p_rho), rho, oo, n=max_order).removeO()
 603                lead = radsimp(simplify(powdenest(s.as_leading_term(rho), force=True)))
 604                power = lead.as_powers_dict().get(rho, None)
 605                coeff = lead / rho**power if power is not None else 0
 606                print("lead =", lead)
 607                print("power =", power)
 608                print("coeff =", coeff)
 609                order = validate_order(power, coeff, [x, y], tol)
 610                if order is not None:
 611                    return order
 612            except Exception as e:
 613                print(f"⚠️ polar expansion failed: {e}")
 614    
 615            try:
 616                print("2D symbol_order - method 2")
 617                z = symbols('z', real=True, positive=True)
 618                xi_eta = {xi: (1/z) * cos(theta), eta: (1/z) * sin(theta)}
 619                p_rho = preprocess_sqrt(self.symbol.subs(xi_eta), 1/z)
 620                s = series(simplify(p_rho), z, 0, n=max_order).removeO()
 621                lead = radsimp(simplify(powdenest(s.as_leading_term(z), force=True)))
 622                power = lead.as_powers_dict().get(z, None)
 623                coeff = lead / z**power if power is not None else 0
 624                print("lead =", lead)
 625                print("power =", power)
 626                print("coeff =", coeff)
 627                order = validate_order(power, coeff, [x, y], tol)
 628                if order is not None:
 629                    return -order
 630            except Exception as e:
 631                print(f"⚠️ fallback z (2D) failed: {e}")
 632            return None
 633    
 634        else:
 635            raise NotImplementedError("Only 1D and 2D supported.")
 636
 637    
 638    def asymptotic_expansion(self, order=3):
 639        """
 640        Compute the asymptotic expansion of the symbol as |ξ| → ∞ (high-frequency regime).
 641    
 642        This method expands the pseudo-differential symbol in inverse powers of the 
 643        frequency variable(s), either in 1D or 2D. It handles both polynomial and 
 644        exponential symbols by performing a series expansion in 1/|ξ| up to the specified order.
 645    
 646        The expansion is performed directly in Cartesian coordinates for 1D symbols.
 647        For 2D symbols, the method uses polar coordinates (ρ, θ) to perform the expansion 
 648        at infinity in ρ, then converts the result back to Cartesian coordinates.
 649    
 650        Parameters
 651        ----------
 652        order : int, optional
 653            Maximum order of the asymptotic expansion. Default is 3.
 654    
 655        Returns
 656        -------
 657        sympy.Expr
 658            The asymptotic expansion of the symbol up to the given order, expressed in Cartesian coordinates.
 659            If expansion fails, returns the original unexpanded symbol.
 660    
 661        Notes:
 662        - In 1D: expansion is performed directly in terms of ξ.
 663        - In 2D: the symbol is first rewritten in polar coordinates (ρ,θ), expanded asymptotically 
 664          in ρ → ∞, then converted back to Cartesian coordinates (ξ,η).
 665        - Handles special case when the symbol is an exponential function by expanding its argument.
 666        - Symbolic normalization is applied early (via `simplify`) for 2D expressions to improve convergence.
 667        - Robust to failures: catches exceptions and issues warnings instead of raising errors.
 668        - Final expression is simplified using `powdenest` and `expand` for improved readability.
 669        """
 670        p = self.symbol
 671    
 672        if self.dim == 1:
 673            xi = symbols('xi', real=True, positive=True)
 674    
 675            try:
 676                # Case: exponential function
 677                if p.func == exp and len(p.args) == 1:
 678                    arg = p.args[0]
 679                    arg_series = series(arg, xi, oo, n=order).removeO()
 680                    expanded = series(exp(expand(arg_series)), xi, oo, n=order).removeO()
 681                    return simplify(powdenest(expanded, force=True))
 682                else:
 683                    expanded = series(p, xi, oo, n=order).removeO()
 684                    return simplify(powdenest(expanded, force=True))
 685    
 686            except Exception as e:
 687                print(f"Warning: 1D expansion failed: {e}")
 688                return p
 689    
 690        elif self.dim == 2:
 691            xi, eta = symbols('xi eta', real=True, positive=True)
 692            rho, theta = symbols('rho theta', real=True, positive=True)
 693    
 694            # Normalize before substitution
 695            p = simplify(p)
 696    
 697            # Substitute polar coordinates
 698            p_polar = p.subs({
 699                xi: rho * cos(theta),
 700                eta: rho * sin(theta)
 701            })
 702    
 703            try:
 704                # Handle exponentials
 705                if p_polar.func == exp and len(p_polar.args) == 1:
 706                    arg = p_polar.args[0]
 707                    arg_series = series(arg, rho, oo, n=order).removeO()
 708                    expanded = series(exp(expand(arg_series)), rho, oo, n=order).removeO()
 709                else:
 710                    expanded = series(p_polar, rho, oo, n=order).removeO()
 711    
 712                # Convert back to Cartesian
 713                norm = sqrt(xi**2 + eta**2)
 714                expansion_cart = expanded.subs({
 715                    rho: norm,
 716                    cos(theta): xi / norm,
 717                    sin(theta): eta / norm
 718                })
 719    
 720                # Final simplifications
 721                result = simplify(powdenest(expansion_cart, force=True))
 722                result = expand(result)
 723                return result
 724    
 725            except Exception as e:
 726                print(f"Warning: 2D expansion failed: {e}")
 727                return p  
 728            
 729    def compose_asymptotic(self, other, order=1, mode='kn', sign_convention=None):
 730        """
 731        Compose two pseudo-differential operators using an asymptotic expansion
 732        in the chosen quantization scheme (Kohn–Nirenberg or Weyl).
 733    
 734        Parameters
 735        ----------
 736        other : PseudoDifferentialOperator
 737            The operator to compose with this one.
 738        order : int, default=1
 739            Maximum order of the asymptotic expansion.
 740        mode : {'kn', 'weyl'}, default='kn'
 741            Quantization mode:
 742            - 'kn' : Kohn–Nirenberg quantization (left-quantized)
 743            - 'weyl' : Weyl symmetric quantization
 744        sign_convention : {'standard', 'inverse'}, optional
 745            Controls the phase factor convention for the KN case:
 746            - 'standard' → (i)^(-n), gives [x, ξ] = +i (physics convention)
 747            - 'inverse' → (i)^(+n), gives [x, ξ] = -i (mathematical adjoint convention)
 748            If None, defaults to 'standard'.
 749    
 750        Returns
 751        -------
 752        sympy.Expr
 753            Symbolic expression for the composed symbol up to the given order.
 754    
 755        Notes
 756        -----
 757        - In 1D (Kohn–Nirenberg):
 758            (p ∘ q)(x, ξ) ~ Σₙ (1/n!) (i sgn)^n ∂_ξⁿ p(x, ξ) ∂_xⁿ q(x, ξ)
 759        - In 1D (Weyl):
 760            (p # q)(x, ξ) = exp[(i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q)] p(x, ξ) q(x, ξ)
 761            truncated at given order.
 762    
 763        Examples
 764        --------
 765        X = a*x, Y = b*ξ
 766        X_op.compose_asymptotic(Y_op, order=3, mode='weyl')
 767        """
 768    
 769        from sympy import diff, factorial, simplify, symbols
 770    
 771        assert self.dim == other.dim, "Operator dimensions must match"
 772        p, q = self.symbol, other.symbol
 773    
 774        # Default sign convention
 775        if sign_convention is None:
 776            sign_convention = 'standard'
 777        sign = -1 if sign_convention == 'standard' else +1
 778    
 779        # --- 1D case ---
 780        if self.dim == 1:
 781            x = self.vars_x[0]
 782            xi = symbols('xi', real=True)
 783            result = 0
 784    
 785            if mode == 'kn':  # Kohn–Nirenberg
 786                for n in range(order + 1):
 787                    term = (1 / factorial(n)) * diff(p, xi, n) * diff(q, x, n) * (1j) ** (sign * n)
 788                    result += term
 789    
 790            elif mode == 'weyl':  # Weyl symmetric composition
 791                # Weyl star product: exp((i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q))
 792                result = 0
 793                for n in range(order + 1):
 794                    for k in range(n + 1):
 795                        # k derivatives acting as (∂_ξ^k p)(∂_x^(n−k) q)
 796                        coeff = (1 / (factorial(k) * factorial(n - k))) * ((1j / 2) ** n) * ((-1) ** (n - k))
 797                        term = coeff * diff(p, xi, k, x, n - k, evaluate=True) * diff(q, x, k, xi, n - k, evaluate=True)
 798                        result += term
 799    
 800            else:
 801                raise ValueError("mode must be either 'kn' or 'weyl'")
 802    
 803            return simplify(result)
 804    
 805        # --- 2D case ---
 806        elif self.dim == 2:
 807            x, y = self.vars_x
 808            xi, eta = symbols('xi eta', real=True)
 809            result = 0
 810    
 811            if mode == 'kn':
 812                for n in range(order + 1):
 813                    for i in range(n + 1):
 814                        j = n - i
 815                        term = (1 / (factorial(i) * factorial(j))) * \
 816                               diff(p, xi, i, eta, j) * diff(q, x, i, y, j) * (1j) ** (sign * n)
 817                        result += term
 818    
 819            elif mode == 'weyl':
 820                for n in range(order + 1):
 821                    for i in range(n + 1):
 822                        j = n - i
 823                        coeff = (1 / (factorial(i) * factorial(j))) * ((1j / 2) ** n) * ((-1) ** (n - i))
 824                        term = coeff * diff(p, xi, i, eta, j, x, 0, y, 0) * diff(q, x, i, y, j, xi, 0, eta, 0)
 825                        result += term
 826            else:
 827                raise ValueError("mode must be either 'kn' or 'weyl'")
 828    
 829            return simplify(result)
 830    
 831        else:
 832            raise NotImplementedError("Only 1D and 2D cases are implemented")
 833
 834    def commutator_symbolic(self, other, order=1, mode='kn', sign_convention=None):
 835        """
 836        Compute the symbolic commutator [A, B] = A∘B − B∘A of two pseudo-differential operators
 837        using formal asymptotic expansion of their composition symbols.
 838    
 839        This method computes the asymptotic expansion of the commutator's symbol up to a given 
 840        order, based on the symbolic calculus of pseudo-differential operators in the 
 841        Kohn–Nirenberg quantization. The result is a purely symbolic sympy expression that 
 842        captures the leading-order noncommutativity of the operators.
 843    
 844        Parameters
 845        ----------
 846        other : PseudoDifferentialOperator
 847            The pseudo-differential operator B to commute with this operator A.
 848        order : int, default=1
 849            Maximum order of the asymptotic expansion. 
 850            - order=1 yields the leading term proportional to the Poisson bracket {p, q}.
 851            - Higher orders include correction terms involving higher mixed derivatives.
 852    
 853        Returns
 854        -------
 855        sympy.Expr
 856            Symbolic expression for the asymptotic expansion of the commutator symbol 
 857            σ([A,B]) = σ(A∘B − B∘A).
 858    
 859        """
 860        assert self.dim == other.dim, "Operator dimensions must match"
 861        p, q = self.symbol, other.symbol
 862    
 863        pq = self.compose_asymptotic(other, order=order, mode=mode, sign_convention=sign_convention)
 864        qp = other.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
 865        
 866        comm_symbol = simplify(pq-qp)
 867
 868        return comm_symbol
 869
 870    def right_inverse_asymptotic(self, order=1):
 871        """
 872        Construct a formal right inverse R of the pseudo-differential operator P such that 
 873        the composition P ∘ R equals the identity plus a smoothing operator of order -order.
 874    
 875        This method computes an asymptotic expansion for the right inverse using recursive 
 876        corrections based on derivatives of the symbol p(x, ξ) and lower-order terms of R.
 877    
 878        Parameters
 879        ----------
 880        order : int
 881            Number of terms to include in the asymptotic expansion. Higher values improve 
 882            approximation at the cost of complexity and computational effort.
 883    
 884        Returns
 885        -------
 886        sympy.Expr
 887            The symbolic expression representing the formal right inverse R(x, ξ), which satisfies:
 888            P ∘ R = Id + O(⟨ξ⟩^{-order}), where ⟨ξ⟩ = (1 + |ξ|²)^{1/2}.
 889    
 890        Notes
 891        -----
 892        - In 1D: The recursion involves spatial derivatives of R and derivatives of p with respect to ξ.
 893        - In 2D: The multi-index generalization is used with mixed derivatives in ξ and η.
 894        - The construction relies on the non-vanishing of the principal symbol p to ensure invertibility.
 895        - Each term in the expansion corresponds to higher-order corrections involving commutators 
 896          between the operator P and the current approximation of R.
 897        """
 898        p = self.symbol
 899        if self.dim == 1:
 900            x = self.vars_x[0]
 901            xi = symbols('xi', real=True)
 902            r = 1 / p.subs(xi, xi)  # r0
 903            R = r
 904            for n in range(1, order + 1):
 905                term = 0
 906                for k in range(1, n + 1):
 907                    coeff = (1j)**(-k) / factorial(k)
 908                    inner = diff(p, xi, k) * diff(R, x, k)
 909                    term += coeff * inner
 910                R = R - r * term
 911        elif self.dim == 2:
 912            x, y = self.vars_x
 913            xi, eta = symbols('xi eta', real=True)
 914            r = 1 / p.subs({xi: xi, eta: eta})
 915            R = r
 916            for n in range(1, order + 1):
 917                term = 0
 918                for k1 in range(n + 1):
 919                    for k2 in range(n + 1 - k1):
 920                        if k1 + k2 == 0: continue
 921                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 922                        dp = diff(p, xi, k1, eta, k2)
 923                        dR = diff(R, x, k1, y, k2)
 924                        term += coeff * dp * dR
 925                R = R - r * term
 926        return R
 927
 928    def left_inverse_asymptotic(self, order=1):
 929        """
 930        Construct a formal left inverse L such that the composition L ∘ P equals the identity 
 931        operator up to terms of order ξ^{-order}. This expansion is performed asymptotically 
 932        at infinity in the frequency variable(s).
 933    
 934        The left inverse is built iteratively using symbolic differentiation and the 
 935        method of asymptotic expansions for pseudo-differential operators. It ensures that:
 936        
 937            L(P(x,ξ),x,D) ∘ P(x,D) = Id + smoothing operator of order -order
 938    
 939        Parameters
 940        ----------
 941        order : int, optional
 942            Maximum number of terms in the asymptotic expansion (default is 1). Higher values 
 943            yield more accurate inverses at the cost of increased computational complexity.
 944    
 945        Returns
 946        -------
 947        sympy.Expr
 948            Symbolic expression representing the principal symbol of the formal left inverse 
 949            operator L(x,ξ). This expression depends on spatial variables and frequencies, 
 950            and includes correction terms up to the specified order.
 951    
 952        Notes
 953        -----
 954        - In 1D: Uses recursive application of the Leibniz formula for symbols.
 955        - In 2D: Generalizes to multi-indices for mixed derivatives in (x,y) and (ξ,η).
 956        - Each term involves combinations of derivatives of the original symbol p(x,ξ) and 
 957          previously computed terms of the inverse.
 958        - Coefficients include powers of 1j (i) and factorial normalization for derivative terms.
 959        """
 960        p = self.symbol
 961        if self.dim == 1:
 962            x = self.vars_x[0]
 963            xi = symbols('xi', real=True)
 964            l = 1 / p.subs(xi, xi)
 965            L = l
 966            for n in range(1, order + 1):
 967                term = 0
 968                for k in range(1, n + 1):
 969                    coeff = (1j)**(-k) / factorial(k)
 970                    inner = diff(L, xi, k) * diff(p, x, k)
 971                    term += coeff * inner
 972                L = L - term * l
 973        elif self.dim == 2:
 974            x, y = self.vars_x
 975            xi, eta = symbols('xi eta', real=True)
 976            l = 1 / p.subs({xi: xi, eta: eta})
 977            L = l
 978            for n in range(1, order + 1):
 979                term = 0
 980                for k1 in range(n + 1):
 981                    for k2 in range(n + 1 - k1):
 982                        if k1 + k2 == 0: continue
 983                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 984                        dp = diff(p, x, k1, y, k2)
 985                        dL = diff(L, xi, k1, eta, k2)
 986                        term += coeff * dL * dp
 987                L = L - term * l
 988        return L
 989
 990    def formal_adjoint(self):
 991        """
 992        Compute the formal adjoint symbol P* of the pseudo-differential operator.
 993
 994        The adjoint is defined such that for any test functions u and v,
 995        ⟨P u, v⟩ = ⟨u, P* v⟩ holds in the distributional sense. This is obtained by 
 996        taking the complex conjugate of the symbol and expanding it asymptotically 
 997        at infinity to ensure proper behavior under integration by parts.
 998
 999        Returns
1000        -------
1001        sympy.Expr
1002            The adjoint symbol P*(x, ξ) in 1D or P*(x, y, ξ, η) in 2D.
1003        
1004        Notes:
1005        - In 1D, the expansion is performed in powers of 1/|ξ|.
1006        - In 2D, the expansion is radial in |ξ| = sqrt(ξ² + η²).
1007        - This method ensures symbolic simplifications for readability and efficiency.
1008        """
1009        p = self.symbol
1010        if self.dim == 1:
1011            x, = self.vars_x
1012            xi = symbols('xi', real=True)
1013            p_star = conjugate(p)
1014            p_star = simplify(series(p_star, xi, oo, n=6).removeO())
1015            return p_star
1016        elif self.dim == 2:
1017            x, y = self.vars_x
1018            xi, eta = symbols('xi eta', real=True)
1019            p_star = conjugate(p)
1020            p_star = simplify(series(p_star, sqrt(xi**2 + eta**2), oo, n=6).removeO())
1021            return p_star
1022
1023    def exponential_symbol(self, t=1.0, order=1, mode='kn', sign_convention=None):
1024        """
1025        Compute the symbol of exp(tP) using asymptotic expansion methods.
1026        
1027        This method calculates the exponential of a pseudo-differential operator 
1028        using either a direct power series expansion or a Magnus expansion, 
1029        depending on the structure of the symbol. The result is valid up to 
1030        the specified asymptotic order.
1031        
1032        Parameters
1033        ----------
1034        t : float or sympy.Symbol, default=1.0
1035            Time or evolution parameter. Common uses:
1036            - t = -i*τ for Schrödinger evolution: exp(-iτH)
1037            - t = τ for heat/diffusion: exp(τΔ)
1038            - t for general propagators
1039        order : int, default=3
1040            Maximum order of the asymptotic expansion. Higher orders include 
1041            more composition terms, improving accuracy for small t or when 
1042            non-commutativity effects are significant.
1043        
1044        Returns
1045        -------
1046        sympy.Expr
1047            Symbolic expression for the exponential operator symbol, computed 
1048            as an asymptotic series up to the specified order.
1049        
1050        Notes
1051        -----
1052        - For commutative symbols (e.g., pure multiplication operators), the 
1053          exponential is exact: exp(tP) = exp(t*p(x,ξ)).
1054        
1055        - For general non-commutative operators, the method uses the BCH-type 
1056          expansion via iterated composition:
1057          exp(tP) ~ I + tP + (t²/2!)P∘P + (t³/3!)P∘P∘P + ...
1058          
1059        - Each power P^n is computed via compose_asymptotic, which accounts 
1060          for the non-commutativity through derivative terms.
1061        
1062        - The expansion is valid for |t| small enough or when the symbol has 
1063          appropriate decay/growth properties.
1064        
1065        - In quantum mechanics (Schrödinger): U(t) = exp(-itH/ℏ) represents 
1066          the time evolution operator.
1067        
1068        - In parabolic PDEs (heat equation): exp(tΔ) is the heat kernel.
1069
1070        """
1071        if self.dim == 1:
1072            x = self.vars_x[0]
1073            xi = symbols('xi', real=True)
1074            
1075            # Initialize with identity
1076            result = 1
1077            
1078            # First order term: tP
1079            current_power = self.symbol
1080            result += t * current_power
1081            
1082            # Higher order terms: (t^n/n!) P^n computed via composition
1083            for n in range(2, order + 1):
1084                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1085                # We use a temporary operator for composition
1086                temp_op = PseudoDifferentialOperator(
1087                    current_power, [x], mode='symbol'
1088                )
1089                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1090                
1091                # Add term (t^n/n!) * P^n
1092                coeff = t**n / factorial(n)
1093                result += coeff * current_power
1094            
1095            return simplify(result)
1096        
1097        elif self.dim == 2:
1098            x, y = self.vars_x
1099            xi, eta = symbols('xi eta', real=True)
1100            
1101            # Initialize with identity
1102            result = 1
1103            
1104            # First order term: tP
1105            current_power = self.symbol
1106            result += t * current_power
1107            
1108            # Higher order terms: (t^n/n!) P^n computed via composition
1109            for n in range(2, order + 1):
1110                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1111                temp_op = PseudoDifferentialOperator(
1112                    current_power, [x, y], mode='symbol'
1113                )
1114                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1115                
1116                # Add term (t^n/n!) * P^n
1117                coeff = t**n / factorial(n)
1118                result += coeff * current_power
1119            
1120            return simplify(result)
1121        
1122        else:
1123            raise NotImplementedError("Only 1D and 2D operators are supported")
1124        
1125    def trace_formula(self, volume_element=None, numerical=False, 
1126                      x_bounds=None, xi_bounds=None):
1127        """
1128        Compute the semiclassical trace of the pseudo-differential operator.
1129        
1130        The trace formula relates the quantum trace of an operator to a 
1131        phase-space integral of its symbol, providing a fundamental link 
1132        between classical and quantum mechanics. This implementation supports 
1133        both symbolic and numerical integration.
1134        
1135        Parameters
1136        ----------
1137        volume_element : sympy.Expr, optional
1138            Custom volume element for the phase space integration. If None, 
1139            uses the standard Liouville measure dx dξ/(2π)^d.
1140        numerical : bool, default=False
1141            If True, perform numerical integration over specified bounds.
1142            If False, attempt symbolic integration (may fail for complex symbols).
1143        x_bounds : tuple of tuples, optional
1144            Spatial integration bounds. For 1D: ((x_min, x_max),)
1145            For 2D: ((x_min, x_max), (y_min, y_max))
1146            Required if numerical=True.
1147        xi_bounds : tuple of tuples, optional
1148            Frequency integration bounds. For 1D: ((xi_min, xi_max),)
1149            For 2D: ((xi_min, xi_max), (eta_min, eta_max))
1150            Required if numerical=True.
1151        
1152        Returns
1153        -------
1154        sympy.Expr or float
1155            The trace of the operator. Returns a symbolic expression if 
1156            numerical=False, or a float if numerical=True.
1157        
1158        Notes
1159        -----
1160        - The semiclassical trace formula states:
1161          Tr(P) = (2π)^{-d} ∫∫ p(x,ξ) dx dξ
1162          where d is the spatial dimension and p(x,ξ) is the operator symbol.
1163        
1164        - For 1D: Tr(P) = (1/2π) ∫_{-∞}^{∞} ∫_{-∞}^{∞} p(x,ξ) dx dξ
1165        
1166        - For 2D: Tr(P) = (1/4π²) ∫∫∫∫ p(x,y,ξ,η) dx dy dξ dη
1167        
1168        - This formula is exact for trace-class operators and provides an 
1169          asymptotic approximation for general pseudo-differential operators.
1170        
1171        - Physical interpretation: the trace counts the "number of states" 
1172          weighted by the observable p(x,ξ).
1173        
1174        - For projection operators (χ_Ω with χ² = χ), the trace gives the 
1175          dimension of the range, related to the phase space volume of Ω.
1176        
1177        - The factor (2π)^{-d} comes from the quantum normalization of 
1178          coherent states / Weyl quantization.
1179        """
1180        from sympy import integrate, simplify, lambdify
1181        from scipy.integrate import dblquad, nquad
1182        
1183        p = self.symbol
1184        
1185        if numerical:
1186            if x_bounds is None or xi_bounds is None:
1187                raise ValueError(
1188                    "x_bounds and xi_bounds must be provided for numerical integration"
1189                )
1190        
1191        if self.dim == 1:
1192            x, = self.vars_x
1193            xi = symbols('xi', real=True)
1194            
1195            if volume_element is None:
1196                volume_element = 1 / (2 * pi)
1197            
1198            if numerical:
1199                # Numerical integration
1200                p_func = lambdify((x, xi), p, 'numpy')
1201                (x_min, x_max), = x_bounds
1202                (xi_min, xi_max), = xi_bounds
1203                
1204                def integrand(xi_val, x_val):
1205                    return p_func(x_val, xi_val)
1206                
1207                result, error = dblquad(
1208                    integrand,
1209                    x_min, x_max,
1210                    lambda x: xi_min, lambda x: xi_max
1211                )
1212                
1213                result *= float(volume_element)
1214                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1215                return result
1216            
1217            else:
1218                # Symbolic integration
1219                integrand = p * volume_element
1220                
1221                try:
1222                    # Try to integrate over xi first, then x
1223                    integral_xi = integrate(integrand, (xi, -oo, oo))
1224                    integral_x = integrate(integral_xi, (x, -oo, oo))
1225                    return simplify(integral_x)
1226                except:
1227                    print("Warning: Symbolic integration failed. Try numerical=True")
1228                    return integrate(integrand, (xi, -oo, oo), (x, -oo, oo))
1229        
1230        elif self.dim == 2:
1231            x, y = self.vars_x
1232            xi, eta = symbols('xi eta', real=True)
1233            
1234            if volume_element is None:
1235                volume_element = 1 / (4 * pi**2)
1236            
1237            if numerical:
1238                # Numerical integration in 4D
1239                p_func = lambdify((x, y, xi, eta), p, 'numpy')
1240                (x_min, x_max), (y_min, y_max) = x_bounds
1241                (xi_min, xi_max), (eta_min, eta_max) = xi_bounds
1242                
1243                def integrand(eta_val, xi_val, y_val, x_val):
1244                    return p_func(x_val, y_val, xi_val, eta_val)
1245                
1246                result, error = nquad(
1247                    integrand,
1248                    [
1249                        [eta_min, eta_max],
1250                        [xi_min, xi_max],
1251                        [y_min, y_max],
1252                        [x_min, x_max]
1253                    ]
1254                )
1255                
1256                result *= float(volume_element)
1257                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1258                return result
1259            
1260            else:
1261                # Symbolic integration
1262                integrand = p * volume_element
1263                
1264                try:
1265                    # Integrate in order: eta, xi, y, x
1266                    integral_eta = integrate(integrand, (eta, -oo, oo))
1267                    integral_xi = integrate(integral_eta, (xi, -oo, oo))
1268                    integral_y = integrate(integral_xi, (y, -oo, oo))
1269                    integral_x = integrate(integral_y, (x, -oo, oo))
1270                    return simplify(integral_x)
1271                except:
1272                    print("Warning: Symbolic integration failed. Try numerical=True")
1273                    return integrate(
1274                        integrand,
1275                        (eta, -oo, oo), (xi, -oo, oo),
1276                        (y, -oo, oo), (x, -oo, oo)
1277                    )
1278        
1279        else:
1280            raise NotImplementedError("Only 1D and 2D operators are supported")
1281
1282    def symplectic_flow(self):
1283        """
1284        Compute the Hamiltonian vector field associated with the principal symbol.
1285
1286        This method derives the canonical equations of motion for the phase space variables 
1287        (x, ξ) in 1D or (x, y, ξ, η) in 2D, based on the Hamiltonian formalism. These describe 
1288        how position and frequency variables evolve under the flow generated by the symbol.
1289
1290        Returns
1291        -------
1292        dict
1293            A dictionary containing the components of the Hamiltonian vector field:
1294            - In 1D: keys are 'dx/dt' and 'dxi/dt', corresponding to dx/dt = ∂p/∂ξ and dξ/dt = -∂p/∂x.
1295            - In 2D: keys are 'dx/dt', 'dy/dt', 'dxi/dt', and 'deta/dt', with similar definitions:
1296              dx/dt = ∂p/∂ξ, dy/dt = ∂p/∂η, dξ/dt = -∂p/∂x, dη/dt = -∂p/∂y.
1297
1298        Notes
1299        -----
1300        - The Hamiltonian here is the principal symbol p(x, ξ) itself.
1301        - This flow preserves the symplectic structure of phase space.
1302        """
1303        if self.dim == 1:
1304            x,  = self.vars_x
1305            xi = symbols('xi', real=True)
1306            return {
1307                'dx/dt': diff(self.symbol, xi),
1308                'dxi/dt': -diff(self.symbol, x)
1309            }
1310        elif self.dim == 2:
1311            x, y = self.vars_x
1312            xi, eta = symbols('xi eta', real=True)
1313            return {
1314                'dx/dt': diff(self.symbol, xi),
1315                'dy/dt': diff(self.symbol, eta),
1316                'dxi/dt': -diff(self.symbol, x),
1317                'deta/dt': -diff(self.symbol, y)
1318            }
1319
1320    def is_elliptic_numerically(self, x_grid, xi_grid, threshold=1e-8):
1321        """
1322        Check if the pseudo-differential symbol p(x, ξ) is elliptic over a given grid.
1323    
1324        A symbol is considered elliptic if its magnitude |p(x, ξ)| remains bounded away from zero 
1325        across all points in the spatial-frequency domain. This method evaluates the symbol on a 
1326        grid of spatial and frequency coordinates and checks whether its minimum absolute value 
1327        exceeds a specified threshold.
1328    
1329        Resampling is applied to large grids to prevent excessive memory usage, particularly in 2D.
1330    
1331        Parameters
1332        ----------
1333        x_grid : ndarray
1334            Spatial grid: either a 1D array (x) or a tuple of two 1D arrays (x, y).
1335        xi_grid : ndarray
1336            Frequency grid: either a 1D array (ξ) or a tuple of two 1D arrays (ξ, η).
1337        threshold : float, optional
1338            Minimum acceptable value for |p(x, ξ)|. If the smallest evaluated symbol value falls below this,
1339            the symbol is not considered elliptic.
1340    
1341        Returns
1342        -------
1343        bool
1344            True if the symbol is elliptic on the resampled grid, False otherwise.
1345        """
1346        RESAMPLE_SIZE = 32  # Reduced size to prevent memory explosion
1347        
1348        if self.dim == 1:
1349            x_vals = x_grid
1350            xi_vals = xi_grid
1351            # Resampling if necessary
1352            if len(x_vals) > RESAMPLE_SIZE:
1353                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1354            if len(xi_vals) > RESAMPLE_SIZE:
1355                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1356        
1357            X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1358            symbol_vals = self.p_func(X, XI)
1359        
1360        elif self.dim == 2:
1361            x_vals, y_vals = x_grid
1362            xi_vals, eta_vals = xi_grid
1363        
1364            # Spatial resampling
1365            if len(x_vals) > RESAMPLE_SIZE:
1366                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1367            if len(y_vals) > RESAMPLE_SIZE:
1368                y_vals = np.linspace(y_vals.min(), y_vals.max(), RESAMPLE_SIZE)
1369        
1370            # Frequency resampling
1371            if len(xi_vals) > RESAMPLE_SIZE:
1372                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1373            if len(eta_vals) > RESAMPLE_SIZE:
1374                eta_vals = np.linspace(eta_vals.min(), eta_vals.max(), RESAMPLE_SIZE)
1375        
1376            X, Y, XI, ETA = np.meshgrid(x_vals, y_vals, xi_vals, eta_vals, indexing='ij')
1377            symbol_vals = self.p_func(X, Y, XI, ETA)
1378        
1379        min_abs_val = np.min(np.abs(symbol_vals))
1380        return min_abs_val > threshold
1381
1382
1383    def is_self_adjoint(self, tol=1e-10):
1384        """
1385        Check whether the pseudo-differential operator is formally self-adjoint (Hermitian).
1386
1387        A self-adjoint operator satisfies P = P*, where P* is the formal adjoint of P.
1388        This property is essential for ensuring real-valued eigenvalues and stable evolution 
1389        in quantum mechanics and symmetric wave propagation.
1390
1391        Parameters
1392        ----------
1393        tol : float
1394            Tolerance for symbolic comparison between P and P*. Small numerical differences 
1395            below this threshold are considered equal.
1396
1397        Returns
1398        -------
1399        bool
1400            True if the symbol p(x, ξ) equals its formal adjoint p*(x, ξ) within the given tolerance,
1401            indicating that the operator is self-adjoint.
1402
1403        Notes:
1404        - The formal adjoint is computed via conjugation and asymptotic expansion at infinity in ξ.
1405        - Symbolic simplification is used to verify equality, ensuring robustness against superficial 
1406          expression differences.
1407        """
1408        p = self.symbol
1409        p_star = self.formal_adjoint()
1410        return simplify(p - p_star).equals(0)
1411
1412    def visualize_fiber(self, x_grid, xi_grid, x0=0.0, y0=0.0):
1413        """
1414        Plot the cotangent fiber structure at a fixed spatial point (x₀[, y₀]).
1415    
1416        This visualization shows how the symbol p(x, ξ) behaves on the cotangent fiber 
1417        above a fixed spatial point. In microlocal analysis, this provides insight into 
1418        the frequency content of the operator at that location.
1419    
1420        Parameters
1421        ----------
1422        x_grid : ndarray
1423            Spatial grid values (1D) for evaluation in 1D case.
1424        xi_grid : ndarray
1425            Frequency grid values (1D) for evaluation in both 1D and 2D cases.
1426        x0 : float, optional
1427            Fixed x-coordinate of the base point in space (1D or 2D).
1428        y0 : float, optional
1429            Fixed y-coordinate of the base point in space (2D only).
1430    
1431        Notes
1432        -----
1433        - In 1D: Displays |p(x, ξ)| over the (x, ξ) phase plane near the fixed point.
1434        - In 2D: Fixes (x₀, y₀) and evaluates p(x₀, y₀, ξ, η), showing the fiber over that point.
1435        - The color map represents the magnitude of the symbol, highlighting regions where it vanishes or becomes singular.
1436    
1437        Raises
1438        ------
1439        NotImplementedError
1440            If called in 2D with missing or improperly formatted grids.
1441        """
1442        if self.dim == 1:
1443            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1444            symbol_vals = self.p_func(X, XI)
1445            plt.contourf(X, XI, np.abs(symbol_vals), levels=50, cmap='viridis')
1446            plt.colorbar(label='|Symbol|')
1447            plt.xlabel('x (position)')
1448            plt.ylabel('ξ (frequency)')
1449            plt.title('Cotangent Fiber Structure')
1450            plt.show()
1451        elif self.dim == 2:
1452            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, xi_grid)
1453            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1454            plt.contourf(xi_grid, xi_grid, np.abs(symbol_vals), levels=50, cmap='viridis')
1455            plt.colorbar(label='|Symbol|')
1456            plt.xlabel('ξ')
1457            plt.ylabel('η')
1458            plt.title(f'Cotangent Fiber at x={x0}, y={y0}')
1459            plt.show()
1460
1461    def visualize_symbol_amplitude(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
1462        """
1463        Display the modulus |p(x, ξ)| or |p(x, y, ξ₀, η₀)| as a color map.
1464    
1465        This method visualizes the amplitude of the pseudodifferential operator's symbol 
1466        in either 1D or 2D spatial configuration. In 2D, the frequency variables are fixed 
1467        to specified values (ξ₀, η₀) for visualization purposes.
1468    
1469        Parameters
1470        ----------
1471        x_grid, y_grid : ndarray
1472            Spatial grids over which to evaluate the symbol. y_grid is optional and used only in 2D.
1473        xi_grid, eta_grid : ndarray
1474            Frequency grids. In 2D, these define the domain over which the symbol is evaluated,
1475            but the visualization fixes ξ = ξ₀ and η = η₀.
1476        xi0, eta0 : float, optional
1477            Fixed frequency values for slicing in 2D visualization. Defaults to zero.
1478    
1479        Notes
1480        -----
1481        - In 1D: Visualizes |p(x, ξ)| over the (x, ξ) grid.
1482        - In 2D: Visualizes |p(x, y, ξ₀, η₀)| at fixed frequencies ξ₀ and η₀.
1483        - The color intensity represents the magnitude of the symbol, highlighting regions where the symbol is large or small.
1484        """
1485        if self.dim == 1:
1486            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1487            symbol_vals = self.p_func(X, XI) 
1488            plt.pcolormesh(X, XI, np.abs(symbol_vals), shading='auto')
1489            plt.colorbar(label='|Symbol|')
1490            plt.xlabel('x')
1491            plt.ylabel('ξ')
1492            plt.title('Symbol Amplitude |p(x, ξ)|')
1493            plt.show()
1494        elif self.dim == 2:
1495            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
1496            XI = np.full_like(X, xi0)
1497            ETA = np.full_like(Y, eta0)
1498            symbol_vals = self.p_func(X, Y, XI, ETA)
1499            plt.pcolormesh(X, Y, np.abs(symbol_vals), shading='auto')
1500            plt.colorbar(label='|Symbol|')
1501            plt.xlabel('x')
1502            plt.ylabel('y')
1503            plt.title(f'Symbol Amplitude at ξ={xi0}, η={eta0}')
1504            plt.show()
1505
1506    def visualize_phase(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
1507        """
1508        Plot the phase (argument) of the pseudodifferential operator's symbol p(x, ξ) or p(x, y, ξ, η).
1509
1510        This visualization helps in understanding the oscillatory behavior and regularity 
1511        properties of the operator in phase space. The phase is displayed modulo 2π using 
1512        a cyclic colormap ('twilight') to emphasize its periodic nature.
1513
1514        Parameters
1515        ----------
1516        x_grid : ndarray
1517            1D array of spatial coordinates (x).
1518        xi_grid : ndarray
1519            1D array of frequency coordinates (ξ).
1520        y_grid : ndarray, optional
1521            2D spatial grid for y-coordinate (in 2D problems). Default is None.
1522        eta_grid : ndarray, optional
1523            2D frequency grid for η (in 2D problems). Not used directly but kept for API consistency.
1524        xi0 : float, optional
1525            Fixed value of ξ for slicing in 2D visualization. Default is 0.0.
1526        eta0 : float, optional
1527            Fixed value of η for slicing in 2D visualization. Default is 0.0.
1528
1529        Notes:
1530        - In 1D: Displays arg(p(x, ξ)) over the (x, ξ) phase plane.
1531        - In 2D: Displays arg(p(x, y, ξ₀, η₀)) for fixed frequency values (ξ₀, η₀).
1532        - Uses plt.pcolormesh with 'twilight' colormap to represent angles from -π to π.
1533
1534        Raises:
1535        - NotImplementedError: If the spatial dimension is not 1D or 2D.
1536        """
1537        if self.dim == 1:
1538            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1539            symbol_vals = self.p_func(X, XI) 
1540            plt.pcolormesh(X, XI, np.angle(symbol_vals), shading='auto', cmap='twilight')
1541            plt.colorbar(label='arg(Symbol) [rad]')
1542            plt.xlabel('x')
1543            plt.ylabel('ξ')
1544            plt.title('Phase Portrait (arg p(x, ξ))')
1545            plt.show()
1546        elif self.dim == 2:
1547            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
1548            XI = np.full_like(X, xi0)
1549            ETA = np.full_like(Y, eta0)
1550            symbol_vals = self.p_func(X, Y, XI, ETA)
1551            plt.pcolormesh(X, Y, np.angle(symbol_vals), shading='auto', cmap='twilight')
1552            plt.colorbar(label='arg(Symbol) [rad]')
1553            plt.xlabel('x')
1554            plt.ylabel('y')
1555            plt.title(f'Phase Portrait at ξ={xi0}, η={eta0}')
1556            plt.show()
1557            
1558    def visualize_characteristic_set(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0, levels=[1e-1]):
1559        """
1560        Visualize the characteristic set of the pseudo-differential symbol, defined as the approximate zero set p(x, ξ) ≈ 0.
1561    
1562        In microlocal analysis, the characteristic set is the locus of points in phase space (x, ξ) where the symbol p(x, ξ) vanishes,
1563        playing a key role in understanding propagation of singularities.
1564    
1565        Parameters
1566        ----------
1567        x_grid : ndarray
1568            Spatial grid values (1D array) for plotting in 1D or evaluation point in 2D.
1569        xi_grid : ndarray
1570            Frequency variable grid values (1D array) used to construct the frequency domain.
1571        x0 : float, optional
1572            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific x position.
1573        y0 : float, optional
1574            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific y position.
1575    
1576        Notes
1577        -----
1578        - For 1D, this method plots the contour of |p(x, ξ)| = ε with ε = 1e-5 over the (x, ξ) plane.
1579        - For 2D, it evaluates the symbol at fixed (x₀, y₀) and plots the characteristic set in the (ξ, η) frequency plane.
1580        - This visualization helps identify directions of degeneracy or hypoellipticity of the operator.
1581    
1582        Raises
1583        ------
1584        NotImplementedError
1585            If called on a solver with dimensionality other than 1D or 2D.
1586    
1587        Displays
1588        ------
1589        A matplotlib contour plot showing either:
1590            - The characteristic curve in the (x, ξ) phase plane (1D),
1591            - The characteristic surface slice in the (ξ, η) frequency plane at (x₀, y₀) (2D).
1592        """
1593        if self.dim == 1:
1594            x_grid = np.asarray(x_grid)
1595            xi_grid = np.asarray(xi_grid)
1596            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1597            symbol_vals = self.p_func(X, XI) 
1598            plt.contour(X, XI, np.abs(symbol_vals), levels=levels, colors='red')
1599            plt.xlabel('x')
1600            plt.ylabel('ξ')
1601            plt.title('Characteristic Set (p(x, ξ) ≈ 0)')
1602            plt.grid(True)
1603            plt.show()
1604        elif self.dim == 2:
1605            if eta_grid is None:
1606                raise ValueError("eta_grid must be provided for 2D visualization.")
1607            xi_grid = np.asarray(xi_grid)
1608            eta_grid = np.asarray(eta_grid)
1609            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
1610            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1611            plt.contour(xi_grid, eta_grid, np.abs(symbol_vals), levels=levels, colors='red')
1612            plt.xlabel('ξ')
1613            plt.ylabel('η')
1614            plt.title(f'Characteristic Set at x={x0}, y={y0}')
1615            plt.grid(True)
1616            plt.show()
1617        else:
1618            raise NotImplementedError("Only 1D/2D characteristic sets supported.")
1619
1620    def visualize_characteristic_gradient(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0):
1621        """
1622        Visualize the norm of the gradient of the symbol in phase space.
1623        
1624        This method computes the magnitude of the gradient |∇p| of a pseudo-differential 
1625        symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. The resulting colormap reveals 
1626        regions where the symbol varies rapidly or remains nearly stationary, 
1627        which is particularly useful for analyzing characteristic sets.
1628        
1629        Parameters
1630        ----------
1631        x_grid : numpy.ndarray
1632            1D array of spatial coordinates for the x-direction.
1633        xi_grid : numpy.ndarray
1634            1D array of frequency coordinates (ξ).
1635        y_grid : numpy.ndarray, optional
1636            1D array of spatial coordinates for the y-direction (used in 2D mode). Default is None.
1637        eta_grid : numpy.ndarray, optional
1638            1D array of frequency coordinates (η) for the 2D case. Default is None.
1639        x0 : float, optional
1640            Fixed x-coordinate for evaluating the symbol in 2D. Default is 0.0.
1641        y0 : float, optional
1642            Fixed y-coordinate for evaluating the symbol in 2D. Default is 0.0.
1643        
1644        Returns
1645        -------
1646        None
1647            Displays a 2D colormap of |∇p| over the relevant phase-space domain.
1648        
1649        Notes
1650        -----
1651        - In 1D, the full gradient ∇p = (∂ₓp, ∂ξp) is computed over the (x, ξ) grid.
1652        - In 2D, the gradient ∇p = (∂ξp, ∂ηp) is computed at a fixed spatial point (x₀, y₀) over the (ξ, η) grid.
1653        - Numerical differentiation is performed using `np.gradient`.
1654        - High values of |∇p| indicate rapid variation of the symbol, while low values typically suggest characteristic regions.
1655        """
1656        if self.dim == 1:
1657            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1658            symbol_vals = self.p_func(X, XI)
1659            grad_x = np.gradient(symbol_vals, axis=0)
1660            grad_xi = np.gradient(symbol_vals, axis=1)
1661            grad_norm = np.sqrt(grad_x**2 + grad_xi**2)
1662            plt.pcolormesh(X, XI, grad_norm, cmap='inferno', shading='auto')
1663            plt.colorbar(label='|∇p|')
1664            plt.xlabel('x')
1665            plt.ylabel('ξ')
1666            plt.title('Gradient Norm (High Near Zeros)')
1667            plt.grid(True)
1668            plt.show()
1669        elif self.dim == 2:
1670            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
1671            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1672            grad_xi = np.gradient(symbol_vals, axis=0)
1673            grad_eta = np.gradient(symbol_vals, axis=1)
1674            grad_norm = np.sqrt(np.abs(grad_xi)**2 + np.abs(grad_eta)**2)
1675            plt.pcolormesh(xi_grid, eta_grid, grad_norm, cmap='inferno', shading='auto')
1676            plt.colorbar(label='|∇p|')
1677            plt.xlabel('ξ')
1678            plt.ylabel('η')
1679            plt.title(f'Gradient Norm at x={x0}, y={y0}')
1680            plt.grid(True)
1681            plt.show()
1682
1683    def plot_hamiltonian_flow(self, x0=0.0, xi0=5.0, y0=0.0, eta0=0.0, tmax=1.0, n_steps=100, show_field=True):
1684        """
1685        Integrate and plot the Hamiltonian trajectories of the symbol in phase space.
1686
1687        This method numerically integrates the Hamiltonian vector field derived from 
1688        the operator's symbol to visualize how singularities propagate under the flow. 
1689        It supports both 1D and 2D problems.
1690
1691        Parameters
1692        ----------
1693        x0, xi0 : float
1694            Initial position and frequency (momentum) in 1D.
1695        y0, eta0 : float, optional
1696            Initial position and frequency in 2D; defaults to zero.
1697        tmax : float
1698            Final integration time for the ODE solver.
1699        n_steps : int
1700            Number of time steps used in the integration.
1701
1702        Notes
1703        -----
1704        - The Hamiltonian vector field is obtained from the symplectic flow of the symbol.
1705        - If the field is complex-valued, only its real part is used for integration.
1706        - In 1D, the trajectory is plotted in (x, ξ) phase space.
1707        - In 2D, the spatial trajectory (x(t), y(t)) is shown along with instantaneous 
1708          momentum vectors (ξ(t), η(t)) using a quiver plot.
1709
1710        Raises
1711        ------
1712        NotImplementedError
1713            If the spatial dimension is not 1D or 2D.
1714
1715        Displays
1716        --------
1717        matplotlib plot
1718            Phase space trajectory(ies) showing the evolution of position and momentum 
1719            under the Hamiltonian dynamics.
1720        """
1721        def make_real(expr):
1722            from sympy import re, simplify
1723            expr = expr.doit(deep=True)
1724            return simplify(re(expr))
1725    
1726        H = self.symplectic_flow()
1727    
1728        if any(im(H[k]) != 0 for k in H):
1729            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
1730    
1731        if self.dim == 1:
1732            x, = self.vars_x
1733            xi = symbols('xi', real=True)
1734    
1735            dxdt_expr = make_real(H['dx/dt'])
1736            dxidt_expr = make_real(H['dxi/dt'])
1737    
1738            dxdt = lambdify((x, xi), dxdt_expr, 'numpy')
1739            dxidt = lambdify((x, xi), dxidt_expr, 'numpy')
1740    
1741            def hamilton(t, Y):
1742                x, xi = Y
1743                return [dxdt(x, xi), dxidt(x, xi)]
1744    
1745            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0], t_eval=np.linspace(0, tmax, n_steps))
1746
1747            if sol.status != 0:
1748                print(f"⚠️ Integration warning: {sol.message}")
1749            
1750            n_points = sol.y.shape[1]
1751            if n_points < n_steps:
1752                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
1753                n_steps = n_points
1754
1755            x_vals, xi_vals = sol.y
1756    
1757            plt.plot(x_vals, xi_vals)
1758            plt.xlabel("x")
1759            plt.ylabel("ξ")
1760            plt.title("Hamiltonian Flow in Phase Space (1D)")
1761            plt.grid(True)
1762            plt.show()
1763    
1764        elif self.dim == 2:
1765            x, y = self.vars_x
1766            xi, eta = symbols('xi eta', real=True)
1767    
1768            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
1769            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
1770            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
1771            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
1772    
1773            def hamilton(t, Y):
1774                x, y, xi, eta = Y
1775                return [
1776                    dxdt(x, y, xi, eta),
1777                    dydt(x, y, xi, eta),
1778                    dxidt(x, y, xi, eta),
1779                    detadt(x, y, xi, eta)
1780                ]
1781    
1782            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0], t_eval=np.linspace(0, tmax, n_steps))
1783
1784            if sol.status != 0:
1785                print(f"⚠️ Integration warning: {sol.message}")
1786            
1787            n_points = sol.y.shape[1]
1788            if n_points < n_steps:
1789                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
1790                n_steps = n_points
1791
1792            x_vals, y_vals, xi_vals, eta_vals = sol.y
1793    
1794            plt.plot(x_vals, y_vals, label='Position')
1795            plt.quiver(x_vals, y_vals, xi_vals, eta_vals, scale=20, width=0.003, alpha=0.5, color='r')
1796            
1797            # Vector field of the flow (optional)
1798            if show_field:
1799                X, Y = np.meshgrid(np.linspace(min(x_vals), max(x_vals), 20),
1800                                   np.linspace(min(y_vals), max(y_vals), 20))
1801                XI, ETA = xi0 * np.ones_like(X), eta0 * np.ones_like(Y)
1802                U = dxdt(X, Y, XI, ETA)
1803                V = dydt(X, Y, XI, ETA)
1804                plt.quiver(X, Y, U, V, color='gray', alpha=0.2, scale=30, width=0.002)
1805
1806            plt.xlabel("x")
1807            plt.ylabel("y")
1808            plt.title("Hamiltonian Flow in Phase Space (2D)")
1809            plt.legend()
1810            plt.grid(True)
1811            plt.axis('equal')
1812            plt.show()
1813
1814    def plot_symplectic_vector_field(self, xlim=(-2, 2), klim=(-5, 5), density=30):
1815        """
1816        Visualize the symplectic vector field (Hamiltonian vector field) associated with the operator's symbol.
1817
1818        The plotted vector field corresponds to (∂_ξ p, -∂_x p), where p(x, ξ) is the principal symbol 
1819        of the pseudo-differential operator. This field governs the bicharacteristic flow in phase space.
1820
1821        Parameters
1822        ----------
1823        xlim : tuple of float
1824            Range for spatial variable x, as (x_min, x_max).
1825        klim : tuple of float
1826            Range for frequency variable ξ, as (ξ_min, ξ_max).
1827        density : int
1828            Number of grid points per axis for the visualization grid.
1829
1830        Raises
1831        ------
1832        NotImplementedError
1833            If called on a 2D operator (currently only 1D implementation available).
1834
1835        Notes
1836        -----
1837        - Only supports one-dimensional operators.
1838        - Uses symbolic differentiation to compute ∂_ξ p and ∂_x p.
1839        - Numerical evaluation is done via lambdify with NumPy backend.
1840        - Visualization uses matplotlib quiver plot to show vector directions.
1841        """
1842        x_vals = np.linspace(*xlim, density)
1843        xi_vals = np.linspace(*klim, density)
1844        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1845
1846        if self.dim != 1:
1847            raise NotImplementedError("Only 1D version implemented.")
1848
1849        x, = self.vars_x
1850        xi = symbols('xi', real=True)
1851        H = self.symplectic_flow()
1852        dxdt = lambdify((x, xi), simplify(H['dx/dt']), 'numpy')
1853        dxidt = lambdify((x, xi), simplify(H['dxi/dt']), 'numpy')
1854
1855        U = dxdt(X, XI)
1856        V = dxidt(X, XI)
1857
1858        plt.quiver(X, XI, U, V, scale=10, width=0.005)
1859        plt.xlabel('x')
1860        plt.ylabel(r'$\xi$')
1861        plt.title("Symplectic Vector Field (1D)")
1862        plt.grid(True)
1863        plt.show()
1864
1865    def visualize_micro_support(self, xlim=(-2, 2), klim=(-10, 10), threshold=1e-3, density=300):
1866        """
1867        Visualize the micro-support of the operator by plotting the inverse of the symbol magnitude 1 / |p(x, ξ)|.
1868    
1869        The micro-support provides insight into the singularities of a pseudo-differential operator 
1870        in phase space (x, ξ). Regions where |p(x, ξ)| is small correspond to large values in 1/|p(x, ξ)|,
1871        highlighting areas of significant operator influence or singularity.
1872    
1873        Parameters
1874        ----------
1875        xlim : tuple
1876            Spatial domain limits (x_min, x_max).
1877        klim : tuple
1878            Frequency domain limits (ξ_min, ξ_max).
1879        threshold : float
1880            Threshold below which |p(x, ξ)| is considered effectively zero; used for numerical stability.
1881        density : int
1882            Number of grid points along each axis for visualization resolution.
1883    
1884        Raises
1885        ------
1886        NotImplementedError
1887            If called on a solver with dimension greater than 1 (only 1D visualization is supported).
1888    
1889        Notes
1890        -----
1891        - This method evaluates the symbol p(x, ξ) over a grid and plots its reciprocal to emphasize 
1892          regions where the symbol is near zero.
1893        - A small constant (1e-10) is added to the denominator to avoid division by zero.
1894        - The resulting plot helps identify characteristic sets.
1895        """
1896        if self.dim != 1:
1897            raise NotImplementedError("Only 1D micro-support visualization implemented.")
1898
1899        x_vals = np.linspace(*xlim, density)
1900        xi_vals = np.linspace(*klim, density)
1901        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1902        Z = np.abs(self.p_func(X, XI))
1903
1904        plt.contourf(X, XI, 1 / (Z + 1e-10), levels=100, cmap='inferno')
1905        plt.colorbar(label=r'$1/|p(x,\xi)|$')
1906        plt.xlabel('x')
1907        plt.ylabel(r'$\xi$')
1908        plt.title("Micro-Support Estimate (1/|Symbol|)")
1909        plt.show()
1910
1911    def group_velocity_field(self, xlim=(-2, 2), klim=(-10, 10), density=30):
1912        """
1913        Plot the group velocity field ∇_ξ p(x, ξ) for 1D pseudo-differential operators.
1914
1915        The group velocity represents the speed at which waves of different frequencies propagate 
1916        in a dispersive medium. It is defined as the gradient of the symbol p(x, ξ) with respect 
1917        to the frequency variable ξ.
1918
1919        Parameters
1920        ----------
1921        xlim : tuple of float
1922            Spatial domain limits (x-axis).
1923        klim : tuple of float
1924            Frequency domain limits (ξ-axis).
1925        density : int
1926            Number of grid points per axis used for visualization.
1927
1928        Raises
1929        ------
1930        NotImplementedError
1931            If called on a 2D operator, since this visualization is only implemented for 1D.
1932
1933        Notes
1934        -----
1935        - This method visualizes the vector field (∂p/∂ξ) in phase space.
1936        - Used for analyzing wave propagation properties and dispersion relations.
1937        - Requires symbolic expression self.expr depending on x and ξ.
1938        """
1939        if self.dim != 1:
1940            raise NotImplementedError("Only 1D group velocity visualization implemented.")
1941
1942        x, = self.vars_x
1943        xi = symbols('xi', real=True)
1944        dp_dxi = diff(self.symbol, xi)
1945        grad_func = lambdify((x, xi), dp_dxi, 'numpy')
1946
1947        x_vals = np.linspace(*xlim, density)
1948        xi_vals = np.linspace(*klim, density)
1949        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1950        V = grad_func(X, XI)
1951
1952        plt.quiver(X, XI, np.ones_like(V), V, scale=10, width=0.004)
1953        plt.xlabel('x')
1954        plt.ylabel(r'$\xi$')
1955        plt.title("Group Velocity Field (1D)")
1956        plt.grid(True)
1957        plt.show()
1958
1959    def animate_singularity(self, xi0=5.0, eta0=0.0, x0=0.0, y0=0.0,
1960                            tmax=4.0, n_frames=100, projection=None):
1961        """
1962        Animate the propagation of a singularity under the Hamiltonian flow.
1963
1964        This method visualizes how a singularity (x₀, y₀, ξ₀, η₀) evolves in phase space 
1965        according to the Hamiltonian dynamics induced by the principal symbol of the operator.
1966        The animation integrates the Hamiltonian equations of motion and supports various projections:
1967        position (x-y), frequency (ξ-η), or mixed phase space coordinates.
1968
1969        Parameters
1970        ----------
1971        xi0, eta0 : float
1972            Initial frequency components (ξ₀, η₀).
1973        x0, y0 : float
1974            Initial spatial coordinates (x₀, y₀).
1975        tmax : float
1976            Total time of integration (final animation time).
1977        n_frames : int
1978            Number of frames in the resulting animation.
1979        projection : str or None
1980            Type of projection to display:
1981                - 'position' : x vs y (or x alone in 1D)
1982                - 'frequency': ξ vs η (or ξ alone in 1D)
1983                - 'phase'    : mixed coordinates like x vs ξ or x vs η
1984                If None, defaults to 'phase' in 1D and 'position' in 2D.
1985
1986        Returns
1987        -------
1988        matplotlib.animation.FuncAnimation
1989            Animation object that can be displayed interactively in Jupyter notebooks or saved as a video.
1990
1991        Notes
1992        -----
1993        - In 1D, only one spatial and one frequency variable are used.
1994        - Complex-valued Hamiltonian fields are truncated to their real parts for integration.
1995        - Trajectories are shown with both instantaneous position (dot) and full path (dashed line).
1996        """
1997        rc('animation', html='jshtml')
1998    
1999        def make_real(expr):
2000            from sympy import re, simplify
2001            expr = expr.doit(deep=True)
2002            return simplify(re(expr))
2003  
2004        H = self.symplectic_flow()
2005
2006        H = {k: v.doit(deep=True) for k, v in H.items()}
2007
2008        print("H = ", H)
2009    
2010        if any(im(H[k]) != 0 for k in H):
2011            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
2012    
2013        if self.dim == 1:
2014            x, = self.vars_x
2015            xi = symbols('xi', real=True)
2016    
2017            dxdt = lambdify((x, xi), make_real(H['dx/dt']), 'numpy')
2018            dxidt = lambdify((x, xi), make_real(H['dxi/dt']), 'numpy')
2019    
2020            def hamilton(t, Y):
2021                x, xi = Y
2022                return [dxdt(x, xi), dxidt(x, xi)]
2023    
2024            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0],
2025                            t_eval=np.linspace(0, tmax, n_frames))
2026            
2027            if sol.status != 0:
2028                print(f"⚠️ Integration warning: {sol.message}")
2029            
2030            n_points = sol.y.shape[1]
2031            if n_points < n_frames:
2032                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2033                n_frames = n_points
2034
2035            x_vals, xi_vals = sol.y
2036    
2037            if projection is None:
2038                projection = 'phase'
2039    
2040            fig, ax = plt.subplots()
2041            point, = ax.plot([], [], 'ro')
2042            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2043    
2044            if projection == 'phase':
2045                ax.set_xlabel('x')
2046                ax.set_ylabel(r'$\xi$')
2047                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2048                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2049    
2050                def update(i):
2051                    point.set_data([x_vals[i]], [xi_vals[i]])
2052                    traj.set_data(x_vals[:i+1], xi_vals[:i+1])
2053                    return point, traj
2054    
2055            elif projection == 'position':
2056                ax.set_xlabel('x')
2057                ax.set_ylabel('x')
2058                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2059                ax.set_ylim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2060    
2061                def update(i):
2062                    point.set_data([x_vals[i]], [x_vals[i]])
2063                    traj.set_data(x_vals[:i+1], x_vals[:i+1])
2064                    return point, traj
2065    
2066            elif projection == 'frequency':
2067                ax.set_xlabel(r'$\xi$')
2068                ax.set_ylabel(r'$\xi$')
2069                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2070                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2071    
2072                def update(i):
2073                    point.set_data([xi_vals[i]], [xi_vals[i]])
2074                    traj.set_data(xi_vals[:i+1], xi_vals[:i+1])
2075                    return point, traj
2076    
2077            else:
2078                raise ValueError("Invalid projection mode")
2079    
2080            ax.set_title(f"1D Singularity Flow ({projection})")
2081            ax.grid(True)
2082            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2083            plt.close(fig)
2084            return ani
2085    
2086        elif self.dim == 2:
2087            x, y = self.vars_x
2088            xi, eta = symbols('xi eta', real=True)
2089    
2090            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
2091            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
2092            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
2093            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
2094    
2095            def hamilton(t, Y):
2096                x, y, xi, eta = Y
2097                return [
2098                    dxdt(x, y, xi, eta),
2099                    dydt(x, y, xi, eta),
2100                    dxidt(x, y, xi, eta),
2101                    detadt(x, y, xi, eta)
2102                ]
2103    
2104            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0],
2105                            t_eval=np.linspace(0, tmax, n_frames))
2106
2107            if sol.status != 0:
2108                print(f"⚠️ Integration warning: {sol.message}")
2109            
2110            n_points = sol.y.shape[1]
2111            if n_points < n_frames:
2112                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2113                n_frames = n_points
2114                
2115            x_vals, y_vals, xi_vals, eta_vals = sol.y
2116    
2117            if projection is None:
2118                projection = 'position'
2119    
2120            fig, ax = plt.subplots()
2121            point, = ax.plot([], [], 'ro')
2122            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2123    
2124            if projection == 'position':
2125                ax.set_xlabel('x')
2126                ax.set_ylabel('y')
2127                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2128                ax.set_ylim(np.min(y_vals) - 1, np.max(y_vals) + 1)
2129    
2130                def update(i):
2131                    point.set_data([x_vals[i]], [y_vals[i]])
2132                    traj.set_data(x_vals[:i+1], y_vals[:i+1])
2133                    return point, traj
2134    
2135            elif projection == 'frequency':
2136                ax.set_xlabel(r'$\xi$')
2137                ax.set_ylabel(r'$\eta$')
2138                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2139                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2140    
2141                def update(i):
2142                    point.set_data([xi_vals[i]], [eta_vals[i]])
2143                    traj.set_data(xi_vals[:i+1], eta_vals[:i+1])
2144                    return point, traj
2145    
2146            elif projection == 'phase':
2147                ax.set_xlabel('x')
2148                ax.set_ylabel(r'$\eta$')
2149                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2150                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2151    
2152                def update(i):
2153                    point.set_data([x_vals[i]], [eta_vals[i]])
2154                    traj.set_data(x_vals[:i+1], eta_vals[:i+1])
2155                    return point, traj
2156    
2157            else:
2158                raise ValueError("Invalid projection mode")
2159    
2160            ax.set_title(f"2D Singularity Flow ({projection})")
2161            ax.grid(True)
2162            ax.axis('equal')
2163            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2164            plt.close(fig)
2165            return ani
2166
2167    def interactive_symbol_analysis(pseudo_op,
2168                                    xlim=(-2, 2), ylim=(-2, 2),
2169                                    xi_range=(0.1, 5), eta_range=(-5, 5),
2170                                    density=100):
2171        """
2172        Launch an interactive dashboard for symbol exploration using ipywidgets.
2173    
2174        This function provides a user-friendly interface to visualize various aspects of the pseudo-differential operator's symbol.
2175        It supports multiple visualization modes in both 1D and 2D, including group velocity fields, micro-support estimates,
2176        symplectic vector fields, symbol amplitude/phase, cotangent fiber structure, characteristic sets and Hamiltonian flows.
2177    
2178        Parameters
2179        ----------
2180        pseudo_op : PseudoDifferentialOperator
2181            The pseudo-differential operator whose symbol is to be analyzed interactively.
2182        xlim, ylim : tuple of float
2183            Spatial domain limits along x and y axes respectively.
2184        xi_range, eta_range : tuple
2185            Frequency domain limits along ξ and η axes respectively.
2186        density : int
2187            Number of points per axis used to construct the evaluation grid. Controls resolution.
2188    
2189        Notes
2190        -----
2191        - In 1D mode, sliders control the fixed frequency (ξ₀) and spatial position (x₀).
2192        - In 2D mode, additional sliders control the second frequency component (η₀) and second spatial coordinate (y₀).
2193        - Visualization updates dynamically as parameters are adjusted via sliders or dropdown menus.
2194        - Supported visualization modes:
2195            'Symbol Amplitude'           : |p(x,ξ)| or |p(x,y,ξ,η)|
2196            'Symbol Phase'               : arg(p(x,ξ)) or similar in 2D
2197            'Micro-Support (1/|p|)'      : Reciprocal of symbol magnitude
2198            'Cotangent Fiber'            : Structure of symbol over frequency space at fixed x
2199            'Characteristic Set'         : Zero set approximation {p ≈ 0}
2200            'Characteristic Gradient'    : |∇p(x, ξ)| or |∇p(x₀, y₀, ξ, η)|
2201            'Group Velocity Field'       : ∇_ξ p(x,ξ) or ∇_{ξ,η} p(x,y,ξ,η)
2202            'Symplectic Vector Field'    : (∇_ξ p, -∇_x p) or similar in 2D
2203            'Hamiltonian Flow'           : Trajectories generated by the Hamiltonian vector field
2204    
2205        Raises
2206        ------
2207        NotImplementedError
2208            If the spatial dimension is not 1D or 2D.
2209    
2210        Prints
2211        ------
2212        Interactive matplotlib figures with dynamic updates based on widget inputs.
2213        """
2214        dim = pseudo_op.dim
2215        expr = pseudo_op.expr
2216        vars_x = pseudo_op.vars_x
2217    
2218        mode_selector_1D = Dropdown(
2219            options=[
2220                'Symbol Amplitude',
2221                'Symbol Phase',
2222                'Micro-Support (1/|p|)',
2223                'Cotangent Fiber',
2224                'Characteristic Set',
2225                'Characteristic Gradient',
2226                'Group Velocity Field',
2227                'Symplectic Vector Field',
2228                'Hamiltonian Flow',
2229            ],
2230            value='Symbol Amplitude',
2231            description='Mode:'
2232        )
2233
2234        mode_selector_2D = Dropdown(
2235            options=[
2236                'Symbol Amplitude',
2237                'Symbol Phase',
2238                'Micro-Support (1/|p|)',
2239                'Cotangent Fiber',
2240                'Characteristic Set',
2241                'Characteristic Gradient',
2242                'Symplectic Vector Field',
2243                'Hamiltonian Flow',
2244            ],
2245            value='Symbol Amplitude',
2246            description='Mode:'
2247        )
2248    
2249        x_vals = np.linspace(*xlim, density)
2250        if dim == 2:
2251            y_vals = np.linspace(*ylim, density)
2252    
2253        if dim == 1:
2254            x, = vars_x
2255            xi = symbols('xi', real=True)
2256            grad_func = lambdify((x, xi), diff(expr, xi), 'numpy')
2257            symplectic_func = lambdify((x, xi), [diff(expr, xi), -diff(expr, x)], 'numpy')
2258            symbol_func = lambdify((x, xi), expr, 'numpy')
2259
2260            xi_slider = FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2261            x_slider = FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2262    
2263            def plot_1d(mode, xi0, x0):
2264                X = x_vals[:, None]
2265    
2266                if mode == 'Group Velocity Field':
2267                    V = grad_func(X, xi0)
2268                    plt.quiver(X, V, np.ones_like(V), V, scale=10, width=0.004)
2269                    plt.xlabel('x')
2270                    plt.title(f'Group Velocity Field at ξ={xi0:.2f}')
2271    
2272                elif mode == 'Micro-Support (1/|p|)':
2273                    Z = 1 / (np.abs(symbol_func(X, xi0)) + 1e-10)
2274                    plt.plot(x_vals, Z)
2275                    plt.xlabel('x')
2276                    plt.title(f'Micro-Support (1/|p|) at ξ={xi0:.2f}')
2277    
2278                elif mode == 'Symplectic Vector Field':
2279                    U, V = symplectic_func(X, xi0)
2280                    plt.quiver(X, V, U, V, scale=10, width=0.004)
2281                    plt.xlabel('x')
2282                    plt.title(f'Symplectic Field at ξ={xi0:.2f}')
2283    
2284                elif mode == 'Symbol Amplitude':
2285                    Z = np.abs(symbol_func(X, xi0))
2286                    plt.plot(x_vals, Z)
2287                    plt.xlabel('x')
2288                    plt.title(f'Symbol Amplitude |p(x,ξ)| at ξ={xi0:.2f}')
2289    
2290                elif mode == 'Symbol Phase':
2291                    Z = np.angle(symbol_func(X, xi0))
2292                    plt.plot(x_vals, Z)
2293                    plt.xlabel('x')
2294                    plt.title(f'Symbol Phase arg(p(x,ξ)) at ξ={xi0:.2f}')
2295    
2296                elif mode == 'Cotangent Fiber':
2297                    pseudo_op.visualize_fiber(x_vals, np.linspace(*xi_range, density), x0=x0)
2298    
2299                elif mode == 'Characteristic Set':
2300                    pseudo_op.visualize_characteristic_set(x_vals, np.linspace(*xi_range, density), x0=x0)
2301    
2302                elif mode == 'Characteristic Gradient':
2303                    pseudo_op.visualize_characteristic_gradient(x_vals, np.linspace(*xi_range, density), x0=x0)
2304    
2305                elif mode == 'Hamiltonian Flow':
2306                    pseudo_op.plot_hamiltonian_flow(x0=x0, xi0=xi0)
2307    
2308            # --- Dynamic container for sliders ---
2309            controls_box = VBox([mode_selector_1D, xi_slider, x_slider])
2310            # --- Function to adjust visible sliders based on mode ---
2311            def update_controls(change):
2312                mode = change['new']
2313                # modes that depend only on xi and eta
2314                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)',
2315                            'Group Velocity Field', 'Symplectic Vector Field']:
2316                    controls_box.children = [mode_selector_1D, xi_slider]
2317                # modes that require xi and x
2318                elif mode in ['Hamiltonian Flow']:
2319                    controls_box.children = [mode_selector_1D, xi_slider, x_slider]
2320                # modes that require nothing
2321                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2322                    controls_box.children = [mode_selector_1D]
2323            mode_selector_1D.observe(update_controls, names='value')
2324            update_controls({'new': mode_selector_1D.value}) 
2325            # --- Interactive binding ---
2326            out = interactive_output(plot_1d, {'mode': mode_selector_1D, 'xi0': xi_slider, 'x0': x_slider})
2327            display(VBox([controls_box, out]))
2328
2329        elif dim == 2:
2330            x, y = vars_x
2331            xi, eta = symbols('xi eta', real=True)
2332            symplectic_func = lambdify((x, y, xi, eta), [diff(expr, xi), diff(expr, eta)], 'numpy')
2333            symbol_func = lambdify((x, y, xi, eta), expr, 'numpy')
2334
2335            xi_slider=FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2336            eta_slider=FloatSlider(min=eta_range[0], max=eta_range[1], step=0.1, value=1.0, description='η₀')
2337            x_slider=FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2338            y_slider=FloatSlider(min=ylim[0], max=ylim[1], step=0.1, value=0.0, description='y₀')
2339    
2340            def plot_2d(mode, xi0, eta0, x0, y0):
2341                X, Y = np.meshgrid(x_vals, y_vals, indexing='ij')
2342    
2343                if mode == 'Micro-Support (1/|p|)':
2344                    Z = 1 / (np.abs(symbol_func(X, Y, xi0, eta0)) + 1e-10)
2345                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='inferno')
2346                    plt.colorbar(label='1/|p|')
2347                    plt.xlabel('x')
2348                    plt.ylabel('y')
2349                    plt.title(f'Micro-Support at ξ={xi0:.2f}, η={eta0:.2f}')
2350    
2351                elif mode == 'Symplectic Vector Field':
2352                    U, V = symplectic_func(X, Y, xi0, eta0)
2353                    plt.quiver(X, Y, U, V, scale=10, width=0.004)
2354                    plt.xlabel('x')
2355                    plt.ylabel('y')
2356                    plt.title(f'Symplectic Field at ξ={xi0:.2f}, η={eta0:.2f}')
2357    
2358                elif mode == 'Symbol Amplitude':
2359                    Z = np.abs(symbol_func(X, Y, xi0, eta0))
2360                    plt.pcolormesh(X, Y, Z, shading='auto')
2361                    plt.colorbar(label='|p(x,y,ξ,η)|')
2362                    plt.xlabel('x')
2363                    plt.ylabel('y')
2364                    plt.title(f'Symbol Amplitude at ξ={xi0:.2f}, η={eta0:.2f}')
2365    
2366                elif mode == 'Symbol Phase':
2367                    Z = np.angle(symbol_func(X, Y, xi0, eta0))
2368                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='twilight')
2369                    plt.colorbar(label='arg(p)')
2370                    plt.xlabel('x')
2371                    plt.ylabel('y')
2372                    plt.title(f'Symbol Phase at ξ={xi0:.2f}, η={eta0:.2f}')
2373    
2374                elif mode == 'Cotangent Fiber':
2375                    pseudo_op.visualize_fiber(np.linspace(*xi_range, density), np.linspace(*eta_range, density),
2376                                              x0=x0, y0=y0)
2377    
2378                elif mode == 'Characteristic Set':
2379                    pseudo_op.visualize_characteristic_set(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2380                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2381    
2382                elif mode == 'Characteristic Gradient':
2383                    pseudo_op.visualize_characteristic_gradient(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2384                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2385    
2386                elif mode == 'Hamiltonian Flow':
2387                    pseudo_op.plot_hamiltonian_flow(x0=x0, y0=y0, xi0=xi0, eta0=eta0)
2388                    
2389            # --- Dynamic container for sliders ---
2390            controls_box = VBox([mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider])
2391            # --- Function to adjust visible sliders based on mode ---
2392            def update_controls(change):
2393                mode = change['new']
2394                # modes that depend only on xi
2395                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)', 'Symplectic Vector Field']:
2396                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider]
2397                # modes that require xi, eta, x and y
2398                elif mode in ['Hamiltonian Flow']:
2399                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider]
2400                # modes that require x and y
2401                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2402                    controls_box.children = [mode_selector_2D, x_slider, y_slider]
2403            mode_selector_2D.observe(update_controls, names='value')
2404            update_controls({'new': mode_selector_2D.value}) 
2405            # --- Interactive binding ---
2406            out = interactive_output(plot_2d, {'mode': mode_selector_2D, 'xi0': xi_slider, 'eta0': eta_slider, 'x0': x_slider, 'y0': y_slider})
2407            display(VBox([controls_box, out]))